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ABSTRACT 

We present the results of a series of radiation-MHD simulations of a local 
patch of an accretion disk, with fixed vertical gravity profile but with differ- 
ent surface mass densities and a broad range of radiation to gas pressure ratios. 
Each simulation achieves a thermal equilibrium that lasts for many cooling times. 
After averaging over times long compared to a cooling time, we find that the ver- 
tically integrated stress is approximately proportional to the vertically-averaged 
total thermal (gas plus radiation) pressure. We map out — for the first time on 
the basis of explicit physics — the thermal equilibrium relation between stress and 
surface density: the stress decreases (increases) with increasing surface mass den- 
sity when the simulation is radiation (gas) pressure dominated. The dependence 
of stress on surface mass density in the radiation pressure dominated regime sug- 
gests the possibility of a Lightman-Eardley inflow instability, but global simula- 
tions or shearing box simulations with much wider radial boxes will be necessary 
to confirm this and determine its nonlinear behavior. 



Subject headings: accretion, accretion disks — instabilities — MHD — X-rays: 
binaries 
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1. Introduction 

It has long been known that models of optically thick, geometric ally thin accretion disks 
based on the alpha stress prescription of IShakura fc Sunyaevl (119731 ) are subject to thermal 
and inflow ( "visco us") instabilities when the vertically- averaged radiation to gas pressure 
ratio exceeds 3/2 (ILightman fc Eardleylll974 IShibazaki &: Hoshil Il975l : IShakura fc Sunyaev 
19761 ). Such radiation pressure dominated accretion disks are expected to be relevant for 
luminous active galactic nuclei and quasars as well as for thermal states of X-ray binaries. 
However, with one possible exception, there has never been clear observational evidence, or 
even observational motivation, for the existence of these instabilities in these sources. This 
is in marked contrast to the cases of dwarf novae and soft X-ray transients, where thermal 
instabilities in the disk associated with h ydrogen ioniz ation, not radiation pressure, are 
central to explaining the observed outbursts (lLasotall200ll ). The X-ray binary GRS 1915+105 
does exhibit recurrent outb urst behavior that has been modelled as being due to radiation 
pressure driven instabilities (IBelloni et al.lll997l ). but it is far from clear that this is the correct 
explanation. This source is the brightest among Galactic black hole X-ray binaries, and 
spend s considerable time at super-Eddington luminosities (IDone. Wardzihski. fc Gierlihski 
20041 ). Other black hole X-ray binaries commonly reach high enough Eddington ratios for 
instabilities to exist according to standard accretion disk theory, but do not exhibit similar 
variability. 

It is now widely suspected that accretion stresses in black hole accretion disks are due 
to turbulence related to the nonlinear growth of the magnetorotational instability (MRI) 
( IBalbus fc Hawleylll998l ). It is computationally feasible to perform thermodynamically con- 
sistent, radiation MHD simulations of this turbulence in local patches of accretion disks. 
These stratified shearing box simulations fully capture grid-scale numerical losses of energy 
as heat and also account for radiative heat losses within the flux-limited diffusion approx- 
imation (IHirose. Krolik. fc Stond 120061 ) . Such simulations have now been performed for a 
broad range of radiation to gas pressure ratios, and in each case an approximate thermal 
equil i brium has been established lasting for many cooling time s (IHirose. Krolik. fc Stone 
20061 : iKrolik. Hirose. & Blaesl l2007i : IHirose. Krolik. fc Blaesl l2009h . No sign of the radiation 
pressure thermal instability is present, even at radia tion to gas pressure ratios we ll above 
the instabilit y threshold of the standard alpha model (IHirose. Krolik. fc Blaesl 120091 ; see also 
Turner! l2004l l 



Thermal stability exists because the stress-pressure relation assumed by the standard 
alpha-model is only established on time scales longer than the thermal timeQ The causal 



In a prescient comment, ILightman fc Eardlevl (|1974l ) suggested that the alpha prescription might only 
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direction of the relation is from stress to pressure, not from pressure to stress. Turbu- 
lence is chaotic and results in a highly fluctuating dissipation rate. It is that dissipation 
that ultimately changes the pressure, but that pressure response is only established after a 
thermal time. An upward fluctuation in pressure does not result in an upward stress re- 
sponse on this time scale as the causal direction is the other way round. Hence there is no 
positive feedback loop on the thermal time scale that would result in a thermal runaway 
rtHirose. Krolik. & Blaesl l2009h . 



(i) 



This still leaves open the question of the slower inflow instability. Mass and angular mo- 
mentum conservation imply that radial mass transport in an accretion disk with a local tur- 
bulen t stress is governed by the equation (jLightman fc Eardleylll974l ; iLynden-Bell fc Pringk 

d£ _ 1 d 
dt r dr 

where £ is the local surface mass density, r is the radius, £' = dljdr is the radial deriva- 
tive of specific angular momentum, and W T< j, is the vertically integrated turbulent stress. 
For geometrically thin disks, the inflow time is much longer than the thermal time, so 
both thermal and vertical hydrostatic equilibrium should be maintained over the time scales 
associated with mass transport. Assuming the disk is optically thick and cools through 
radiative losses (which is the case for all of our strat i fied shearing box simulations so far ; 
Hirose. Krolik. fc Stone! bood : iKrolik. Hirose. fc Blaesl 120071 : Iflirose. Krolik. fc Blaesl l2009h . 
then radiative, hydrostatic, and thermal equilibrium imply that the vertically-averaged stress 
in a radiation pressure dominated disk is simply given by 

eft 2 

where ft is the angular velocity in t he disk, c is the speed of light, and ft' = dQ/dr is 
the shear (IShakura fc Sunyaevl Il976l ) . The opacity k is generally dominated by electron 
scattering in this regime, and is therefore constant. Hydrostatic equilibrium implies that 
the disk half-thickness is H ~ 2P/(ft 2 £), where P is the midplane pressure. Hence at a 
particular radiation pressure dominated radius, 

W r4> ~ 2Hr r s cx ^. (3) 



A standard alpha disk model with r ri p = aP implies from equation (T5]) that P is independent 
of E, so diWrtf)) I dYj < 0. Equation ([1]) therefore represents a diffusion equation with a nega- 
tive diffusion coefficien t. Unstable growth of surface density enhanc ements and rarefactions 
would therefore result ( Lightman fc Eardley 1974 ; Lightman 1974al lbh . 



be valid on slow time scales, which they identified as being of order the inflow time and longer. 
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Like the fictitious thermal instability, however, the reality of the inflow instabil ity has 



always been questionable. For example, as pointed out by lLightman fc Eardleyl (119741 ) them- 
selves, a stress proportional to the gas pressure alone would produce an inflow stable disk. 
On the other hand, the stratified shearing box simulations appear to be consistent with total 



therm al pressure scaling with stress on supra-thermal time scales (IHirose. Krolik. fc Blaes 



20091 ). a fact that we will demonstrate much more explicitly in this paper. 



A plot of vertically integrated stress versus surface density for a range of ther- 
mal equilibria at a fixed radius within a disk would suggest inflow stability or instability 
depending on whether the slope is positive or negative, respectively. Because thermal equi- 
librium implies that the local radiation flux F emerging from each face of the disk is given 
by F = W r( j ) r\fl'\/2, a plot of F (or effective temperature) versus surface density may be 
used in the same way. Such "S-curve" plots are commonly used to investigate t he hydrogen 



ionization-dr iven disk instabilities in dwarf novae and soft X-ray transients (e.g. ISmak!ll984 



Lasotal 120011 ) . 



Using radiation MHD simulations of stratified shearing boxes, we have now mapped 
out the stress versus surface density thermal equilibrium curve for a wide range of radiation 
to gas pressure ratios at a fixed radius in a disk around a stellar mass black hole. This is 
the first time that this curve has been drawn on the basis of explicit physical mechanisms, 
rather than phenomenological estimates. While there are large fluctuations which produce 
an inherent scatter in the thermal equilibrium curve, the results are consistent with the 
standard alpha disk model. The inflow instability might therefore be present even in MRI 
turbulent disks. 

This paper is organized as follows. In section 2 we provide a brief overview of the 
numerical parameters and properties of the simulations. In section 3 we discuss the stress- 
pressure relation and demonstrate that average total thermal pressure, rather than a pressure 
which singles out gas pressure as being special in some way, is best correlated with average 
stress. In section 4, we summarize the resulting thermal equilibria on a stress-surface density 
diagram, and discuss the possible implications. We summarize our conclusions in section 5. 



Simulations 



The radiation MHD equations for stratified shearin g boxes, and the numerical met hods 
we use to solve them, have be en described in detail by IHirose. Krolik. &: Stone! (120061 ) and 
Hirose. Krolik. &: Blaesl (120091 ). and references therein. Grid scale losses of mechanical and 
magnetic energy are fully captured as heat in the gas, and gas and radiation exchange heat 
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through Planck mean free-free absorption and emission and Compton scattering. Radiation 
transport is treated through flux-limited diffusion. 

All the simulations were run with an angular velocity Q = 190 s^ 1 , corresponding to 
a radius of 30GM/c 2 around a 6.62 M Schwarzschild black hole. Different total surface 
mass densities were chosen for each simulation in order to map out the stress-surface density 
relation. Table [T] summarizes the parameters of each simulation. The x, y, and z axes 
correspond to the radial, azimuthal, and vertical directions, respectively. The simulations 
were initialized in hydrostatic and thermal equilibrium under an assumed initial vertical 
profile of dissipation per unit volume proportional to the density divided by the square root 
of the optical depth measured from the nearest surface. Once the simulation starts, the 
dissipation is thereafter self-consistently determined from the turbulent dynamics. Each 
simulation started with a weak magnetic field consistin g of a twisted azimuthal flux t ube 



located at the center of the box. We refer the reader to iHirose. Krolik. fc Blaed (120091 ) for 
more details. 

The ver tical box height L 7 of the sim ulations was chosen so that two conditions would 



be satisfied ( IHirose. Krolik. fc Blaed l2009i ). First, the total surface mass density changes by 
less than a few percent due to vertical mass loss and mass creation by the density floor of 
the simulation. Second, the MRI is always well-resolved in the midplane regions. These 
conditions were checked a posteriori. The upper and lower photospheres are always within 
the simulation domain, although the optically thin regions do not necessarily have a tremen- 
dous vertical extent. This should not significantly affect the emergent flux, which is what 
really matters. Because the radiation energy density in the optically thin regions is nearly 
independent of height, the radiation energy density at the photospheres should depend only 
weakly on box height. It is this radiation energy density that provides the effective boundary 
condition for the radiation diffusion equation in the optically thick regions where the vast 
majority of the dissipation occurs, so the emergent flux should be well-determined. Nev- 
ertheless, we warn the reader that we have not demonstrated numerical convergence with 
respect to variations in the simulation box dimensions. 

Six of the simulations (0211b, 0519b, 1112a, 1126b, 0520a, and 0320a) were initialized 
as radiation pressure dominated, while two (090304a and 090423a) started out gas pressure 
dominated. (Due to an initialization error, simulations 090304a and 090423a actually had 
an angular velocity 1.6% smaller than the others: 187 rad s _1 rather than 190 rad s _1 . 
We do not believe this significantly affects our conclusions, as the fluctuations in stress and 
radiative cooling are considerably l arger than this.) The r e sults of simulations 1112a and 



1 1 26b were discussed extensively in IHirose. Krolik. fc Blaed ( 120091 ) 



All of the simulations share many of the properties that we discussed in detail in previous 
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pape r s on radiation MHD simulations of stratified shearing boxes (IHirose. Krolik. fc Stone 



2006 



Krolik. Hirose. fc Blaesl 120071 : Blaes. Hirose. fc Krolikl 120071 : IHirose. Krolik. fc Blaes 



20091 ). The subphotospheric regions consist of a magnetorotational turbulent zone in the 
midplane regions which is supported against gravity by gas and radiation pressure gradients. 
Further out, magnetic forces dominate or contribute substantially, and Parker instability 
dynamics, rather than MRI turbulence, appear to control the structure of the outer layers. 

Figure [T] shows the thermal and turbulent energy content in the box as a function of 
time for each of the radiation pressure dominated simulations. The thermal history of the 
more gas-dominated simulations are shown in Figure [2j Defining the instantaneous thermal 
time as the total radiation and gas internal energy divided by the emergent radiative flux 
on both vertical faces of the box, the thermal time averaged over the duration of each of 
the simulations ranges from a minimum of 6 orbits for 090423a to a maximum of 24 orbits 
for 0519b. All the simulations have reached an approximate thermal equilibrium, albeit 
with long time scale fluctuati ons. There is no evidence for the thermal instability predicted 
by classic alpha disk models (IShakura fc Sunyaevlll976l ). in spite of the fact that the time- 
averaged ratio of vertically-averaged radiation pressure to vertically-averaged gas pressure is 
as high as 70 in the case of simulation 0519b. 



3. The Stress-Pressure Relation 



A number of authors have suggested alternative stress prescriptions in which the ac- 
cretion stress is proportional to gas pressure or some combination of gas and radiation 
pressures, rather than total thermal pressure (gas plus radiation), in part to stabilize the 
radiation pressure dominated portion of black hole accretion disks (ISakimoto fc Coroniti 



1981 



Stella &: Rosner 



1984 



Burm 1985: Szuszkiewicz 1990; Merloni Sz Fabian 



2002: 



Merloni 



2003 ) Hsakimoto fc Coronitil ~ 1989 ) have even argued that standard alpha disk models are 



inconsistent in the radiation pressure dominated regime, as magnetic fields that are strong 
enough to provide accretion stress would be too buoyant to be retained by the disk. Their 
argument has two flaws, however. First, they assumed that the magnetic field consisted of 
discrete flux tubes, rather than being more continuously distributed throughout the plasma. 
Second, they were unaware at the time of magnetic field generation by the MRI. Note that 
there is no indication in the energy histories shown in Figured] that the magnetic and turbu- 
lent kinetic energies are limited by the gas internal energy. Indeed, in the two most radiation 
pressure dominated simulations (0211b and 0519b), the magnetic energy almost always ex- 
ceeds the gas internal energy, and even the turbulent kinetic energy can occasionally be 
larger than the gas internal energy. 



- 7- 



Based on the behavior of the stress to thermal pressure ratio within simulation 1112a, 
we argued in the past that a prescription where stress and total thermal pressure are propor- 
tional to each other (on time scales longer than the thermal time) is a superior description 
of the simulation behavior than alternatives where the pressure is taken to be just the 
gas pressure alone or the ge ometric mean of radiation and gas pressure (see Fig. 10 of 



Hirose. Krolik. fc Blaea 120091 ). A similar conclusion can be drawn when data is compared 



across the simulations we are considering here. 

For each of our simulations, we computed the vertically-averaged stress r r( ^ av and the 
box-average of various measures of pressure P av : the sum of the radiation pressure P ra d 
and gas pressure P gas , the geometric mean (P ra dPgas) 1,/2 of these pressures, and just the gas 
pressure. We then computed the time average of the ratio of these spatial averages of stress 
and pressure, a =< T r )av /P av >, and compared them to the time average of the ratio of 
spatially averaged radiation and gas pressure, < P ra d,av/Pgas,av >• These time averages were 
computed over the entire duration of each simulation, excluding the first 10 orbits as the 
MRI was still in its growth phase. 

The results are shown in Figure [3J The black points show the results for P av defined 
as the total pressure, and are clearly closest to being independent of the radiation to gas 
pressure ratio. (A linear fit to the dependence of these values of a on the radiation to 
gas pressure ratio gives a slightly negative slope which is consistent with zero within one 
standard deviation of the slope determination.) A weighted average of these values gives 
a = 0.018 ±0.002. The green curve shows what one would obtain if stress and total pressure 
were proportional with this ratio, but we redefined a in terms of the geometric mean of 
radiation and gas pressure, i.e. 



' rm,av 



P 



a 



p \ 1/2 / p \ -1/2 

r rad,av \ / -'rad^v 

P / \ P 

eras p,\r I \ - 1 c 



av / \ - 1 gas,av / \ * gas,av 



(4) 



The blue curve shows the same thing for the gas pressure stress prescription, 



P f-4 i+ (^-))- (5) 

av / \ \ r gas,av / / 

These curves clearly explain the trends seen in the data for these alternative stress prescrip- 
tions. 

Our results are therefore most consistent with the standard alpha prescription involv- 
ing total thermal p ressure. Recent axisymmetric global radiation MHD simulations by 



Ohsuga et aJj (120091 ) also reach a similar conclusion. We emphasize, however, that such 



a prescription should still be treated with caution. There is no obvious reason that a should 
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be a universal constant, and we do not know what physics is determining its value in our 
simulations. It is possible, for example, that a is a function of fl or M, as we have not 
varied these parameters at all; global simulations lHawley fc Kroliki (120011 . |2002| ) have shown 
that it can be a function of radius when the underlying orbital dynamics change. In fact, 
it is rather surprising that it is as const ant as it is when defined in te rms of total thermal 
pressure. Moreover, as we emphasized in lHirose. Krolik. fc Blaed (120091 ). the fact that stress 
is in any way related to pressure is really because turbulent dissipation heats the plasma. 
Only when averaged over height and averaged over many thermal times does the standard 
alpha prescription provide an adequate description of the fact that pressure is correlated 
with stress. 



4. The Stress-Surface Density Relation 

For each simulation, we computed the vertically integrated stress as a function of time, 
and then time averaged this over the simulation's duration, again ignoring the first ten orbits. 
The results are plotted as a function of surface density in Figure HI Vertical error bars indicate 
the standard deviations of the instantaneous fluctuations in stress about the mean. (Due to 
mass loss and an imposed density floor, the surface density also fluctuates, but by at most 
two percent in all the simulations.) The right hand axis indicates the effective temperature 
of the radiation leaving each vertical face of the box if perfect thermal equilibrium held. We 
have also computed versions of the diagram by time-averaging the radiation flux leaving both 
faces of the box and computing the effective temperature, and time-averaging the volume- 
integrated dissipation rate. All three versions are almost identical, as they must be, given 
the approximate thermal equilibrium that has been established in each simulation. 

The curves in Figure H] show the results predicted by standard alpha disk models with 
stress scaling with total pressure, but with internal structure parameters based on the vertical 
structure observed in the simulations themselves. These parameters are defined in Appendix 
A, and listed in Table[2J The simulation data are clearly consistent with the alpha disk model, 
although there are variations in the internal structure parameters that, when averaged over 
all the simulations, produce an alpha disk curve that misses many of the data points (solid 
curve in Figure H]). Note that the alpha disk model predicts a maximum value of surface 
density S crit above which an alpha disk cannot be in thermal equilibrium. We have not 
attempted to simulate such high surface densities. 

The negative slope of the stress-surface density relation in the radiation pressure domi- 
nated simulation data indicates a negative mass diffusion coefficient, suggesting inflow insta- 
bility if local thermal equilibrium is everywhere maintained on the time scales associated with 



- 9- 



nonlocal, radial mass transport. Because this would be a diffusive instability, a pe rturbation 
with radial wavelength A would have a characteristic growth rate ~ aQ(H/X) 2 (jLightman 
1974bl ; IShakura Sunyaevi Il976l ). hence growing faster at shorte r wavele ngths. However, 
there must be a short wavelength cutoff to this trend. iLightmanl (Il974bl ) argues that this 
cutoff i s of order the vertical scal e height H because of turbulent mixing on such length 
scales. IShakura fc Sunyaevi (119761 ) reach approximately the same conclusion based on the 
fact that the linear thermal and inflow unstable modes of the alpha disk equations become 
degenerate and then stable for radial wavelengths less than of order the disk scale height. 
This latter fact suggests that the short wavelength cutoff to the inflow instability might 
actually be considerably larger than the disk scale height, because the radiation pressure 
driven thermal instability itself does not manifest itself in real MRI turbulence. Insofar as 
thermal effects become increasingly important to the inflow instability at short wavelengths, 
this may make the cutoff wavelength longer. 

Because the radial extent of our simulation domain is so small (significantly less than 
a scale height in all cases - see Table [T|), we do not expect (and do not find) the inflow 
instability to be present in our simulations^] In principle, vertically stratified shearing boxes 
with much wider radial domains could produce the instability, even though the shearing box 
has no net radial accretion. In spite of our use of the term "inflow" instability, the instability 
is really one of negative diffusion, and can grow even in the absence of a net mass accretion 
flow through the box. We can demonstrate this fact explicitly by vertically integrating the 
shearing box equations (see Appendix B). This results in a radial mass diffusion equation, 



~dt 



1 <9 2 



(6) 



identical to the local limit of equation ([I]), where q = — dlnQ/dlnr and \x\ <C r. If 
dWxy/dT, < under conditions of thermal equilibrium, we would then expect inflow in- 
stability. Stratified shearing box simulations sufficiently wide in the radial direction might 
therefore manifest this instability in the radiation pressure dominated regime. To see this 
effect may require radial box widths much larger than the vertical pressure scale height, in 



2 Because the thermal instability of the alpha disk also has a short radial wavelength cutoff, the reader 
might wonder how we can claim with confidence that t he thermal instability does no t exist based only on 
shearing box simulations with narrow radial domains ( Hirose. Krolik. fc Blaes 20091 ). This is because an 



infinite radial wavelength mode (one that does not vary at all in radius) can still be present in such boxes, 
and would still grow if it were truly thermally unstable. Physically, the pure thermal instability of the 
standard alpha disk arises because the vertically integrated heating rate is more temperature sensitive than 
the cooling rate, and radial variations are irrelevant. In contrast, the growth rate of infinite radial wavelength 
modes on the inflow instability branch is zero. 
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order that different radial regions be thermally decoupled and that stress scale with local 
pressure on time scales longer than the local thermal time. In any case, such simulations 
would enable a determination of whether the inflow instability is real, and what its minimum 
and fastest growing wavelengths are. 

If the instability does manifest itself in MRI turbulence, such simulations would also help 
determine its nonlinear outcome. Of course, a number of hydrodynamic simulations have 
been done of the combined thermal /inflow instability of radiation dominated alpha-disks 



e.g. iHonma. Matsumoto. fc Katolll99ll ; ISzuszkiewicall998l ). However, the time evolution of 



these simulations is dominated by the faster thermal instability, which we now know to be 
fictitious. As far as we are aware, the only simulations that have ever been attempt ed of the 



radiation dominated inflow instability on its own were done by iLightmanl (Il974bl ) himself, 
who numerically solved the alpha disk mass diffusion equation, assuming that thermal equi- 
librium is strictly maintained. All his simulations therefore started with a surface density 
less than S crit . The resulting evolution rapidly produced clumping of the surface density up 
to S cr it with optically thin rarefied regions in between. Both of these conditions violated 
the assumptions on which t he mass d i ffusion equation is based, and the numerical calcula- 



tion had to be stopped. As ILightmanl (Il974al ) pointed out, surface densities exceeding E crit 



cannot be in thermal equilibrium within the assumptions of the alpha model, because heat 
generation always exceeds cooling through vertical radiative diffusion. On the other hand, 
if the characteristic radial size of the clumps is as small as the disk scale height, radial heat 
transport is likely to be important in the nonlinear outcome. Radiation MHD simulations 
with radially wide shearing boxes could address all these issues. 

Radiation MHD simulations could also clear up another question clouding prediction 
of the outcome of this putative instability: the negative slope of the stress-surface density 
relation on the radiation dominated branch implies growing surface density fluctuations only 
if local thermal equilibrium with purely vertical heat flow is maintained everywhere while 
the instability tries to develop. It could be that subtleties in the behavior of MRI turbulence 
preclude this from happening, just as the thermodynamics of the turbulence prevented the 
thermal instability from manifesting itself. 

This issue is closely related to the question of exactly what the gas-dominated and 
radiation-dominated branches of equilibria in Figure H] truly represent. If this were a low- 
dimensional dynamical system like the standard alpha-disk, then the fact that both branches 
are thermally stable would imply an unstable equilibrium branch in between. However, a 
stratified shearing box with real turbulence is not a low-dimensional dynamical system — the 
spatial dependence of gas density, pressure, velocity, magnetic field, and radiation pressure 
all influence its evolution, and Figure 0] is probably better viewed as a projection of a very 
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complicated dynamical phase space. In our experience, simulations that are initialized far 
away from the equilibrium branches do not undergo steady heating or cooling, but instead 
wander chaotically due to the fluctuating character of the turbulence. Whether there is a 
true unstable thermal equilibrium branch between the gas and radiation-dominated branches 
would be very well masked by this stochasticity even if the vertically-integrated stress were 
the only significant dynamical variable; the much higher dimensionality of the real system 
makes it essentially impossible to test whether such a branch exists on the basis of simulation 
data. It is conceivable that a local perturbation in surface density would require considerable 
time to reach the thermal equilibrium branch, and would in any case fluctuate about that 
branch. It is an open question whether or not this evolution would be conducive to inflow 
instability. 

It is also noteworthy that hydrostatic and radiative equilibrium necessarily enforce a 



chara cteristic stress on the radiation dominated branch (equation [2] IShakura fc Sunyaev 



19761 ). There is no such constraint on the gas dominated branch, and in fact hydrostatic 
equilibrium is not even needed to derive the relation between stress and surface density in the 
alpha model on this branch (see Appendix A). It is possible that these differences may also 
be relevant to evolution on the inflow time scale in the gas and radiation pressure dominated 
regimes in real MRI turbulence. 



5. Conclusions 

We have completed vertically stratified, local radiation MHD simuations of magnetorota- 
tional turbulence with fixed vertical gravity over a range of surface mass densities. All of the 
simulations reach a thermal equilibrium, but with conti nued long term fluctuations in the in 



tenia! energy content. We have confirmed earlier work (jTurnerll2004l ; iHirose. Krolik. &: Blaes 



20091 ) that the radiation pressure dominated thermal instability predicted by the standard 
alpha disk model does not exist, even though the the box- and time- averaged radiation to 
gas pressure ratios in the new simulations are as high as 70. 

However, we also find that, when averaged over many thermal times, the vertically 
integrated total thermal pressure (i.e., radiation plus gas pressure) is well-correlated with the 
vertically integrated stress. Neither the vertically integrated gas pressure nor the geometric 
mean of gas and radiation pressure exhibit such a good correlation. The same simulation 
data therefore yield a thermal equilibrium relation between surface density and (long) time- 
averaged vertically integrated stress that follows the one predicted by the traditional alpha 
model: (W r( p) oc E _1 . Consequently, if thermal equilibrium is maintained on the long time 
scales associated with radial mass transport by the turbulent stresses, these local results 
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sugges t that the nonlocal clumping of mass associated with the classic iLightman fc Eardley 



( 119741 ) inflow instability might develop. Radiation MHD shearing box simulations with much 
wider boxes in the radial direction, or well-resolved global simulations, will be necessary to 
investigate this possibility. 
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A. The Stress/Surface Density Relation for Alpha-Disks Accounting for 

Internal Structure from the Simulations 

We derive one-zone vertical structure equations for an alpha-disk in terms of parameters 
measured from the simulations as follows. Vertical hydrostatic equilibrium implies that the 
midplane pressure is given by 

P(0) = ^ 2 £# pl , (Al) 

where H p i is a density scale height defined through the first vertical moment of the density 
distribution, 

2 r°° 

H pl = -J^ p{z)zdz. (A2) 
The midplane pressure in the simulations is always dominated by gas and radiation pressure, 

^4 ar <°> 1 + w (A3) 

where T(0) is the midplane temperature and fi = 0.6 atomic mass units is the mean particle 
mass assumed in the simulations. We have eliminated the midplane density by defining the 
zeroth vertical moment of the density distribution, 

X ^ POO 

H ^m=mh p(z)dz - (A4) 



The fundamental assumption of an alpha-disk model is that the vertically integrated 
stress Wrtf) is a times the vertical integral of the thermal pressure P t h- We may therefore 
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write 

W r4> = 2aH P P(0), (A5) 
where the thermal pressure scale height has been defined as 



OG 



H P = —J— / P th (z)dz (A6) 



2P(0) J 

Thermal equilibrium implies that the flux F = crT e 4 ff emerging from each face of the disk is 
given by 

F = - l -W r4> T d ^ = -nW r<t> (A7) 
2 dr 4 

Finally, the fact that most of the accretion power ultimately escapes vertically through 
radiative diffusion in the simulations motivates us to write the emergent flux as 

where the opacity is dominated by electron scattering in our simulations, k ~ 0.33. We 
have introduced a parameter £ which is usually taken to be approximately unity in alpha- 
disk models. As we discuss below, however, our measured values of £ in the simulation are 
substantially greater than unity. This is presumably due to the fact that the dissipation 
profile peaks off the midplane, and that a non-negligible fraction of the accretion power is 
transported away from the midplane by mechanical motions. 

Equations OAip . (1A3[) . (1A5[) . (1A7[) . and f1A8h can be combined to give the relationship 
between emergent flux or vertically integrated stress and surface density. In the gas and 
radiation pressure dominated limits, the result is 



3 5 k ( aUlA 4 ( _ffp X '' 



1/3 

\ — ) S5/3 ' if P gas(0) > Prad(0); 

&e(^)^\ ifP r ad(0)»P gas (0). 



(A9) 



This gives the usual result that the alpha disk is inflow stable if gas pressure dominated, 
but inflow unstable if radiation pressure dominated. Note that the hydrostatic equilibrium 
equation (I All) is not needed to derive the gas pressure dominated relation. 

If neither pressure dominates at the midplane, the following equation can be used to 
derive the relationship between flux or stress and surface density: 

5j p3/4 _ 2 ^5/4 _ 3j p5/4£l/2 = ^ (AlQ) 
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where £ = S/S crit , F = F/F 



crit; 



J crit 



>19 



(9)5 lo n 



OLK 



V\ 4 H pl H pO 
k 



H p 



1/8 



and 



12c 2 t 2 Q (H. 



crit 



I pi 



(AH) 



(A12) 



Apart from the extra dimensionless factors of £ and the pressure and density scale heights, 
equation (lAllO agrees with the expression for the maximum su rface density consistent with 
thermal equilibrium that was first derived by iLightmanl (11974al ) . 

Equations (1A9j) and (lAlOh are what we use to plot the alpha disk predictions in Figure HI 
We measured the dimensionless parameters describing the internal structure of the disk 
as follows. Using the time series of horizontally-averaged data from each simulation, we 
measured a as the ratio of the vertically integrated stress to vertically-averaged pressure, 
and £ in terms of the ratio of emerging flux (averaged over both faces of the disk) to midplane 
radiation energy density We then averaged both of these quantities over time. We measured 
the scale height ratios H p0 /Hp and H p i/Hp from the time and horizontally-averaged vertical 
profiles of density and thermal pressure from each simulation. The results are shown in Table 

El 

The scale height ratio parameters are remarkably constant across all the simulations. 
The numerical values of these parameters are close to what one would obtain from simple 
analytic equilibria in a gravitational field that increases linearly with height. An n = 3 
polytrope (adiabatic, radiation pressure dominated) has H pQ /Hp = 1.125 and H p \jHp ~ 
0.67, while an n = 3/2 polytrope (adiabatic, gas pressure dominated) has H p0 /H P = 1.2 
and H p i/Hp ~ 0.69. Note that both ratios are slightly larger in the gas pressure dominated 
polytrope, in agreement with the trend that we see in the simulations. These simple models 
do not agree exactly with the simulations due in part to nonzero entropy gradients. In the 
simulations, the entropy generally increases away from the midplane, so that the density 
must decrease faster for a given pressure decrease compared to an adiabatic profile. Hence 
the ratio of density scale height to pressure scale height is smaller than that of an adiabatic 
profile, and this is why the ratios we measure in the simulations are slightly smaller than 
the polytropic ratios. 

In contrast to the scale height ratios, the stress parameter a shows a little more scat- 
ter. The parameter £ clearly increases as the simulations become more radiation pressure 
dominated. 
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B. The Radial Mass Diffusion Equation for the Shearing Box 



The radial mass diffusion equation can be derived from the equations of the shearing 
box ( Hawley. Gammie. fc Balbud Il995l ) as follows. Define the azimuthal averaged surface 
density at some radius x and time t as 



E(x,t) = 



1 



x+Ax/2 



L y Ax J x _Ax/2 



dx' 



Ly/2 



dy 



L z /2 



dzp(x',y,z,t), 



(Bl) 



'-Ly/2 J-L z /2 

and the mass-weighted vertical and azimuthal average of some quantity X as 



< X > (x,t) 



LyAxTj J x ^Ax/2 



x+Ax/2 



dx' 



Ly/2 



Ly/2 



dy 



L z /2 



L z /2 



dzp(x',y,z,t)X(x',y,z,t). 



:B2) 



Here we have also performed a radial average over a length scale Aa; of order an assumed 
radial coherence length of the turbulence, which is presumably of order the disk scale height. 
Apart from the use of Cartesian coordinates in a box, these definitions are exactly the same as 
the vertical i ntegrations and averages used to derive the alpha-disk equations from the MHD 
equations by lBalbus fc Papaloizoul (119991 ) . Applying this averaging and vertical integration 
to the mass continuity and y-momentum equations of the shearing box, we obtain 



m + o~x^ <v ^ >) 







and 



9 \ 



d s dW xv 

_ (E <„ t >< % » + ^-» 



-2fiS < V r >, 



(B3) 



(B4) 



where W xy is the vertically-integrated and azimuthally-averaged Maxwell and Reynolds 

stress, 

1 



W 



, y\X, t) 

Lin, 



j-Ly/2 


rLz/2 i 


/ dy 


\ dz 


'-Ly/2 J 


-L z /2 V 



B X By 

47T 



pv x Sv, 



(B5) 



Radial derivatives are defined through differences over the length scale Ax. We have also 
assumed zero mass and Poynting flux through the vertical boundaries, and < 5v y >= 0. 
Equations (1B3I) and (1B4I) can then be combined into a radial mass diffusion equation, 



i d 2 w 



(B6) 



where q = — dlnQ/dlnr is the shear parameter. 
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Fig. 1. — Time history of the box integrated radiation internal energy (red), gas internal 
energy (green), magnetic energy (blue) and turbulent kinetic energy (black) for each of the 
radiation pressure dominated simulations. 
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Fig. 2. — Same as Figure [T] except for simulations 090304a and 090423a, which are not 
radiation pressure dominated. Note the much smaller energies on the vertical axis. 
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Fig. 3. — Measured values of the stress parameter a as a function of the time-averaged ratio 
of the box-averaged radiation pressure to the box-averaged gas pressure. The black points 
define a as the time-averaged ratio of the vertically averaged stress to the box-averaged 
total thermal pressure. The blue and green points define a in the same way except with 
total thermal pressure replaced by gas pressure and the geometric mean of gas and radiation 
pressure, respectively. Both horizontal and vertical error bars indicate one standard deviation 
in the time averages. The horizontal black line indicates the weighted mean of a for the total 
pressure stress prescription. Assuming that stress really scales with total thermal pressure, 
the green and blue curves show how a would then behave under the other stress prescriptions. 
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Fig. 4. — Time averaged, vertically integrated stress as a function of surface mass density 
for each simulation. The right hand axis shows the corresponding effective temperature of 
the radiation leaving each vertical face of the box. The solid curve shows the prediction 
of the vertically integrated alpha disk model with internal structure parameters averaged 
over all the simulations. The dashed line shows the prediction of the radiation pressure 
dominated alpha disk model with internal structure parameters averaged over the radia- 
tion pressure dominated simulations. The dot-dashed line shows the prediction of the gas 
pressure dominated alpha disk model with internal structure parameters measured from the 
gas pressure dominated simulations. (See Appendix A for the equations used to define the 
internal structure parameters.) 
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Table 1. 


Simulation Parameters 




Simulation 


E 




H 


Duration 


Box Dimensions 


Grid Zones 




(g cm" 




(cm) 


(orbits) 


{L x /H x L y /H x L z /H) 


(N x xN y x N z ) 


0211b 


5.43 x 


10 4 


5.83 x 10 6 


264 


0.3375 x 1.35 x 6.3 


48 x 96 x 896 


0519b 


7.48 x 


10 4 


4.37 x 10 6 


403 


0.3375 x 1.35 x 6.3 


48 x 96 x 896 


1112a 


1.06 x 


10 5 


1.46 x 10 6 


610 


0.45 x 1.8 x 8.4 


48 x 96 x 896 


1126b 


1.06 x 


10 5 


1.46 x 10 6 


611 


0.45 x 1.8 x 8.4 


48 x 96 x 896 


0520a 


1.24 x 


10 5 


1.17 x 10 6 


603 


0.54 x 2.16 x 10.08 


48 x 96 x 896 


0320a 


1.52 x 


10 5 


7.28 x 10 5 


426 


0.6 x 2.4 x 11.2 


48 x 96 x 896 


090304a 


5.00 x 


10 4 


3.13 x 10 5 


600 


0.625 x 2.5 x 10 


32 x 64 x 512 


090423a 


2.00 x 


10 4 


2.10 x 10 5 


262 


0.5 x 2.0 x 8.0 


32 x 64 x 512 



Table 2. Internal Structure Parameters Measured from Simulations 



Simulation E a £ H p0 /Hp H p i/Hp 

(g cm" 2 ) 



0211b 


5.43 


x 


10 4 


0.019 ±0.005 


5.4 ± 


1.2 


0.952 


0.629 


0519b 


7.48 


x 


10 4 


0.015 ±0.004 


5.6 ± 


1.3 


0.949 


0.630 


1112a 


1.06 


x 


10 5 


0.023 ±0.007 


4.5 ± 


1.1 


0.934 


0.624 


1126b 


1.06 


x 


10 5 


0.020 ±0.006 


4.7 ± 


1.0 


0.948 


0.632 


0520a 


1.24 


x 


10 5 


0.016 ±0.006 


4.6 ± 


1.0 


0.926 


0.629 


0320a 


1.52 


x 


10 5 


0.015 ±0.004 


4.1 ± 


0.8 


0.990 


0.645 


090304a 


5.00 


x 


10 4 


0.028 ±0.009 


2.8 ± 


0.6 


1.090 


0.680 


090423a 


2.00 


x 


10 4 


0.028 ± 0.007 


2.3 ± 


0.4 


1.127 


0.685 


Mean rad a 








0.017 ±0.002 


4.7 ± 


0.4 


0.950 


0.631 


Mean gas b 








0.028 ±0.006 


2.5 ± 


0.4 


1.109 


0.683 


Mean a n c 








0.018 ±0.002 


3.4 ± 


0.3 


0.990 


0.644 



a Averaged over the six radiation pressure dominated simulations, i.e. excluding 090304a 
and 090423a. 

b Averaged over the the gas pressure dominated simulations 090304a and 090423a. 
c Averaged over all eight simulations. 



